# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating figures and tables in Appendix I: Negative Binomial and Count Outcome Models
# Appendix I Table 27: Neg. Binomial Model: White/Non-White Racial Geographic Boundaries Are Associated With Felony Arrests
# Appendix I Table 28: Binary Outome Model: White/Non-White Racial Geographic Boundaries Are Associated With Misdemeanor Arrests
# Appendix I Table 29: Binary Outome Model: White/Non-White Racial Geographic Boundaries Are Associated With Felony Arrests
# Appendix I Table 30: Neg. Binomial Model: Crime Moderation Analysis for White/Non-White Boundaries and Misdemeanor Arrests
# Appendix I Table 31: Neg. Binomial Model: Crime Moderation Analysis for White/Non-White Boundaries and Felony Arrests
# Appendix I Table 32: Binary Outcome: Crime Moderation Analysis for White/Non-White Boundaries and Misdemeanor Arrests
# Appendix I Table 33: Binary Outcome Model: Crime Moderation Analysis for White/Non-White Boundaries and Felony Arrests
# Appendix I Table 34: Neg. Binomial Model: Influence of Logged Crime on Logged Stops (Standardized), Conditional on Racial Boundary Status.
# Appendix I Table 35: Binary Outome Model: Influence of Logged Crime on Logged Stops (Standardized), Conditional on Racial Boundary Status.
# Appendix I Table 36: Neg. Binomial Model: Influence of Percent White on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix I Table 37: Neg. Binomial Model: Influence of Percent White on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix I Table 38: Binary Outome Model: Influence of Percent White on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix I Table 39: Binary Outome Model: Influence of Percent White on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix I Table 40: Negative Binomial Model: Influence of Percent White on Police Stops by Race of Civilian, Conditional on Racial Boundary Status
# Appendix I Table 41: Binary Outcome Model: Influence of Percent White on Police Stops by Race of Civilian, Conditional on Racial Boundary Status


suppressPackageStartupMessages(
  
  {
    library(AER)
    library(dplyr)
    library(fixest)
    library(lfe)
    library(tidyverse)
    library(ggplot2)
    library(MASS)
    library(sensemakr)
    library(haven)
    library(readstata13)
    library(readxl)
    library(readr)
    library(gridExtra)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggthemes)
    library(meta)
  }
)

# first prepare variables by city
# read in data by city

# Load Atlanta final dataset
load("atl_final.RData")

# log and +1 to variables (pop already logged)
atl_fin = atl_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
atl_fin$larrest_sd <- atl_fin$larrests - (mean(atl_fin$larrests)/sd(atl_fin$larrests))
atl_fin$lmisdemeanors_sd <- atl_fin$lmisdemeanors - (mean(atl_fin$lmisdemeanors)/sd(atl_fin$lmisdemeanors))
atl_fin$lfelonies_sd <- atl_fin$lfelonies - (mean(atl_fin$lfelonies)/sd(atl_fin$lfelonies))
atl_fin$lnonviolent_sd <- atl_fin$lnonviolent - (mean(atl_fin$lnonviolent)/sd(atl_fin$lnonviolent))
atl_fin$lviolent_sd <- atl_fin$lviolent - (mean(atl_fin$lviolent)/sd(atl_fin$lviolent))
atl_fin$lsociety_sd <- atl_fin$lsociety - (mean(atl_fin$lsociety)/sd(atl_fin$lsociety))
atl_fin$lperson_sd <- atl_fin$lperson - (mean(atl_fin$lperson)/sd(atl_fin$lperson))
atl_fin$lproperty_sd <- atl_fin$lproperty - (mean(atl_fin$lproperty)/sd(atl_fin$lproperty))

# rename racial and economic boundary variable
atl_fin$atl_white_blv <- atl_fin$p_race_white_blv
atl_fin$atl_ses_blv <- atl_fin$ses_blv

# Load Austin final dataset
load("aus_final.RData")



# log and +1 to variables (pop already logged)
aus_fin = aus_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_fin$larrest_sd <- aus_fin$larrests - (mean(aus_fin$larrests)/sd(aus_fin$larrests))
aus_fin$lmisdemeanors_sd <- aus_fin$lmisdemeanors - (mean(aus_fin$lmisdemeanors)/sd(aus_fin$lmisdemeanors))
aus_fin$lfelonies_sd <- aus_fin$lfelonies - (mean(aus_fin$lfelonies)/sd(aus_fin$lfelonies))
aus_fin$lnonviolent_sd <- aus_fin$lnonviolent - (mean(aus_fin$lnonviolent)/sd(aus_fin$lnonviolent))
aus_fin$lviolent_sd <- aus_fin$lviolent - (mean(aus_fin$lviolent)/sd(aus_fin$lviolent))
aus_fin$lsociety_sd <- aus_fin$lsociety - (mean(aus_fin$lsociety)/sd(aus_fin$lsociety))
aus_fin$lperson_sd <- aus_fin$lperson - (mean(aus_fin$lperson)/sd(aus_fin$lperson))
aus_fin$lproperty_sd <- aus_fin$lproperty - (mean(aus_fin$lproperty)/sd(aus_fin$lproperty))

# rename racial and economic boundary variable
aus_fin$aus_white_blv <- aus_fin$p_race_white_blv
aus_fin$aus_ses_blv <- aus_fin$ses_blv


## Load Boston data
load("bos_final.RData")


# log and +1 to variables (pop already logged)
bos_fin = bos_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime + 1))



# create scaled DVs (to account for different pop. sizes and density in blocks)
bos_fin$larrests_sd <- bos_fin$larrests - (mean(bos_fin$larrests)/sd(bos_fin$larrests))
bos_fin$lmisdemeanors_sd <- bos_fin$lmisdemeanors - (mean(bos_fin$lmisdemeanors)/sd(bos_fin$lmisdemeanors))
bos_fin$lfelonies_sd <- bos_fin$lfelonies - (mean(bos_fin$lfelonies)/sd(bos_fin$lfelonies))
bos_fin$lnonviolent_sd <- bos_fin$lnonviolent - (mean(bos_fin$lnonviolent)/sd(bos_fin$lnonviolent))
bos_fin$lviolent_sd <- bos_fin$lviolent - (mean(bos_fin$lviolent)/sd(bos_fin$lviolent))
bos_fin$lsociety_sd <- bos_fin$lsociety - (mean(bos_fin$lsociety)/sd(bos_fin$lsociety))
bos_fin$lperson_sd <- bos_fin$lperson - (mean(bos_fin$lperson)/sd(bos_fin$lperson))
bos_fin$lproperty_sd <- bos_fin$lproperty - (mean(bos_fin$lproperty)/sd(bos_fin$lproperty))

# rename racial and economic boundary variable
bos_fin$bos_white_blv <- bos_fin$p_race_white_blv
bos_fin$bos_ses_blv <- bos_fin$ses_blv

## Load Chicago data
load("chi_final.RData")

# log and +1 to variables (pop already logged)
chi_fin = chi_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime+1),
         lviolentcrime = log(violent_crime+1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_fin$larrest_sd <- chi_fin$larrests - (mean(chi_fin$larrests)/sd(chi_fin$larrests))
chi_fin$lmisdemeanors_sd <- chi_fin$lmisdemeanors - (mean(chi_fin$lmisdemeanors)/sd(chi_fin$lmisdemeanors))
chi_fin$lnonviolent_sd <- chi_fin$lnonviolent - (mean(chi_fin$lnonviolent)/sd(chi_fin$lnonviolent))
chi_fin$lsociety_sd <- chi_fin$lsociety - (mean(chi_fin$lsociety)/sd(chi_fin$lsociety))
chi_fin$lfelonies_sd <- chi_fin$lfelonies - (mean(chi_fin$lfelonies)/sd(chi_fin$lfelonies))
chi_fin$lperson_sd <- chi_fin$lperson - (mean(chi_fin$lperson)/sd(chi_fin$lperson))
chi_fin$lproperty_sd <- chi_fin$lproperty - (mean(chi_fin$lproperty)/sd(chi_fin$lproperty))
chi_fin$lviolent_sd <- chi_fin$lviolent - (mean(chi_fin$lviolent)/sd(chi_fin$lviolent))


# rename racial and economic boundary variable
chi_fin$chi_white_blv <- chi_fin$p_race_white_blv
chi_fin$chi_ses_blv <- chi_fin$ses_blv

## Load Louisville data
load("lou_final.RData")


# log and +1 to variables (pop already logged) (no misdemeanors vs felonies for louisville)
lou_fin = lou_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property+1),
         lviolentcrime = log(crime_violent+1),
         lviolent = log(violent_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
lou_fin$larrest_sd <- lou_fin$larrests - (mean(lou_fin$larrests)/sd(lou_fin$larrests))
lou_fin$lnonviolent_sd <- lou_fin$lnonviolent - (mean(lou_fin$lnonviolent)/sd(lou_fin$lnonviolent))
lou_fin$lsociety_sd <- lou_fin$lsociety - (mean(lou_fin$lsociety)/sd(lou_fin$lsociety))
lou_fin$lviolent_sd <- lou_fin$lviolent - (mean(lou_fin$lviolent)/sd(lou_fin$lviolent))
lou_fin$lperson_sd <- lou_fin$lperson - (mean(lou_fin$lperson)/sd(lou_fin$lperson))
lou_fin$lproperty_sd <- lou_fin$lproperty - (mean(lou_fin$lproperty)/sd(lou_fin$lproperty))

# rename racial and economic boundary variable
lou_fin$lou_white_blv <- lou_fin$p_race_white_blv
lou_fin$lou_ses_blv <- lou_fin$ses_blv

## Load Milwaukee data
load("mil_final.RData")


# log and +1 to variables (pop already logged)
mil_fin = mil_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_fin$larrest_sd <- mil_fin$larrests - (mean(mil_fin$larrests)/sd(mil_fin$larrests))
mil_fin$lmisdemeanors_sd <- mil_fin$lmisdemeanors - (mean(mil_fin$lmisdemeanors)/sd(mil_fin$lmisdemeanors))
mil_fin$lfelonies_sd <- mil_fin$lfelonies - (mean(mil_fin$lfelonies)/sd(mil_fin$lfelonies))
mil_fin$lnonviolent_sd <- mil_fin$lnonviolent - (mean(mil_fin$lnonviolent)/sd(mil_fin$lnonviolent))
mil_fin$lviolent_sd <- mil_fin$lviolent - (mean(mil_fin$lviolent)/sd(mil_fin$lviolent))
mil_fin$lsociety_sd <- mil_fin$lsociety - (mean(mil_fin$lsociety)/sd(mil_fin$lsociety))
mil_fin$lperson_sd <- mil_fin$lperson - (mean(mil_fin$lperson)/sd(mil_fin$lperson))
mil_fin$lproperty_sd <- mil_fin$lproperty - (mean(mil_fin$lproperty)/sd(mil_fin$lproperty))

# rename racial and economic boundary variable
mil_fin$mil_white_blv <- mil_fin$p_race_white_blv
mil_fin$mil_ses_blv <- mil_fin$ses_blv


## Load Seattle data
load("sea_final.RData")

# log and +1 to variables (pop already logged)
sea_fin = sea_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
sea_fin$larrest_sd <- sea_fin$larrests - (mean(sea_fin$larrests)/sd(sea_fin$larrests))
sea_fin$lmisdemeanors_sd <- sea_fin$lmisdemeanors - (mean(sea_fin$lmisdemeanors)/sd(sea_fin$lmisdemeanors))
sea_fin$lfelonies_sd <- sea_fin$lfelonies - (mean(sea_fin$lfelonies)/sd(sea_fin$lfelonies))
sea_fin$lnonviolent_sd <- sea_fin$lnonviolent - (mean(sea_fin$lnonviolent)/sd(sea_fin$lnonviolent))
sea_fin$lviolent_sd <- sea_fin$lviolent - (mean(sea_fin$lviolent)/sd(sea_fin$lviolent))
sea_fin$lsociety_sd <- sea_fin$lsociety - (mean(sea_fin$lsociety)/sd(sea_fin$lsociety))
sea_fin$lperson_sd <- sea_fin$lperson - (mean(sea_fin$lperson)/sd(sea_fin$lperson))
sea_fin$lproperty_sd <- sea_fin$lproperty - (mean(sea_fin$lproperty)/sd(sea_fin$lproperty))

# rename racial and economic boundary variable
sea_fin$sea_white_blv <- sea_fin$p_race_white_blv
sea_fin$sea_ses_blv <- sea_fin$ses_blv


#### robustness: using different outcome measures to circumvent problem with log outcome, main results #### 

# extensive margin

atl_fin$mis_arr_bin = ifelse(atl_fin$misdemeanor_arrests > 0, 1, 0)
aus_fin$mis_arr_bin = ifelse(aus_fin$misdemeanor_arrests > 0, 1, 0)
bos_fin$mis_arr_bin = ifelse(bos_fin$misdemeanor_arrests > 0, 1, 0)
chi_fin$mis_arr_bin = ifelse(chi_fin$misdemeanor_arrests > 0, 1, 0)
mil_fin$mis_arr_bin = ifelse(mil_fin$misdemeanor_arrests > 0, 1, 0)
sea_fin$mis_arr_bin = ifelse(sea_fin$misdemeanor_arrests > 0, 1, 0)

atl_fin$fel_arr_bin = ifelse(atl_fin$felony_arrests > 0, 1, 0)
aus_fin$fel_arr_bin = ifelse(aus_fin$felony_arrests > 0, 1, 0)
bos_fin$fel_arr_bin = ifelse(bos_fin$felony_arrests > 0, 1, 0)
chi_fin$fel_arr_bin = ifelse(chi_fin$felony_arrests > 0, 1, 0)
mil_fin$fel_arr_bin = ifelse(mil_fin$felony_arrests > 0, 1, 0)
sea_fin$fel_arr_bin = ifelse(sea_fin$felony_arrests > 0, 1, 0)

# misdemeanor arrests (poisson)
# atlanta: null
# austin: positive
# boston: positive
# chicago: null
# milwaukee: null
# seattle: null

# overdispersion test! 

disptest1 = dispersiontest(glm(misdemeanor_arrests ~ 1,
                               data = subset(atl_fin, pop > 0),
                               family = 'poisson'), trafo = 1)
disptest2 = dispersiontest(glm(misdemeanor_arrests ~ 1,
                               data = subset(aus_fin, pop > 0),
                               family = 'poisson'), trafo = 1)
disptest3 = dispersiontest(glm(misdemeanor_arrests ~ 1,
                               data = subset(bos_fin, pop > 0),
                               family = 'poisson'), trafo = 1)
disptest4 = dispersiontest(glm(misdemeanor_arrests ~ 1,
                               data = subset(chi_fin, pop > 0),
                               family = 'poisson'), trafo = 1)
disptest5 = dispersiontest(glm(misdemeanor_arrests ~ 1,
                               data = subset(mil_fin, pop > 0),
                               family = 'poisson'), trafo = 1)
disptest6 = dispersiontest(glm(misdemeanor_arrests ~ 1,
                               data = subset(sea_fin, pop > 0),
                               family = 'poisson'), trafo = 1)

# felony arrests (nbinom)

negbinm1 = fixest::fenegbin(misdemeanor_arrests ~ atl_white_blv + 
                              atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                              lpropertycrime + lviolentcrime + phys_bound + com_density,
                            data = subset(atl_fin, pop > 0), cluster = ~BG_CODE)

negbinm2 = fixest::fenegbin(misdemeanor_arrests ~ aus_white_blv +
                              aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                              lpropertycrime + lviolentcrime + phys_bound + com_density,
                            data = subset(aus_fin, pop > 0), cluster = ~BG_CODE)

negbinm3 = fixest::fenegbin(misdemeanor_arrests ~ bos_white_blv + 
                              bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                            data = bos_fin, cluster = ~BG_CODE)

negbinm4 = fixest::fenegbin(misdemeanor_arrests ~ chi_white_blv + 
                              chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                              lpropertycrime + lviolentcrime + phys_bound + com_density, 
                            data = subset(chi_fin, pop > 0), cluster = ~BG_CODE)

negbinm5 = fixest::fenegbin(misdemeanor_arrests ~ mil_white_blv + 
                              mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + 
                              pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                            data = subset(mil_fin, pop > 0), cluster = ~BG_CODE)

negbinm6 = fixest::fenegbin(misdemeanor_arrests ~ sea_white_blv + 
                              sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                              lpropertycrime + lviolentcrime + phys_bound + com_density, 
                            data = subset(sea_fin, pop > 0), cluster = ~BG_CODE)

# felony arrests (nbinom)

negbinm7 = 
  fixest::fenegbin(felony_arrests ~ atl_white_blv + atl_ses_blv + 
                     lpop + p_race_white + age_15_35_male + 
                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed +
                     pcol + lpropertycrime + lviolentcrime + phys_bound + com_density,
                   data = subset(atl_fin, pop > 0), cluster = ~BG_CODE)

negbinm8 = 
  fixest::fenegbin(felony_arrests ~ aus_white_blv + aus_ses_blv + 
                     lpop + p_race_white + age_15_35_male + 
                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed + 
                     pcol + lpropertycrime + lviolentcrime + phys_bound + com_density,
                   data = subset(aus_fin, pop > 0), cluster = ~BG_CODE)

negbinm9 = 
  fixest::fenegbin(felony_arrests ~ bos_white_blv + bos_ses_blv + 
                     lpop + p_race_white + age_15_35_male + 
                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                   data = bos_fin, cluster = ~BG_CODE)

negbinm10 = 
  fixest::fenegbin(felony_arrests ~ chi_white_blv + chi_ses_blv + 
                     lpop + p_race_white + age_15_35_male + 
                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed + 
                     pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                   data = subset(chi_fin, pop > 0), cluster = ~BG_CODE)

negbinm11 = 
  fixest::fenegbin(felony_arrests ~ mil_white_blv + mil_ses_blv + 
                     lpop + p_race_white + age_15_35_male + 
                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed + 
                     pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                   data = subset(mil_fin, pop > 0), cluster = ~BG_CODE)

negbinm12 = 
  fixest::fenegbin(felony_arrests ~ sea_white_blv + sea_ses_blv + 
                     lpop + p_race_white + age_15_35_male + 
                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed + 
                     pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                   data = subset(sea_fin, pop > 0), cluster = ~BG_CODE)


# binary model

binmod1 = lm_robust(mis_arr_bin ~ atl_white_blv + 
                      atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(atl_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')

binmod2 = lm_robust(mis_arr_bin ~ aus_white_blv +
                      aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(aus_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')

binmod3 = lm_robust(mis_arr_bin ~ bos_white_blv + 
                      bos_ses_blv + lpop + p_race_black + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                    data = bos_fin, cluster = BG_CODE, se_type = 'stata')

binmod4 = lm_robust(mis_arr_bin ~ chi_white_blv + 
                      chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                      lpropertycrime + lviolentcrime + phys_bound + com_density, 
                    data = subset(chi_fin, pop > 0), cluster = BG_CODE)

binmod5 = lm_robust(mis_arr_bin ~ mil_white_blv + 
                      mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + 
                      pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                    data = subset(mil_fin, pop > 0), cluster = BG_CODE)

binmod6 = lm_robust(mis_arr_bin ~ sea_white_blv + 
                      sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                      lpropertycrime + lviolentcrime + phys_bound + com_density, 
                    data = subset(sea_fin, pop > 0), cluster = BG_CODE)


binmod7 = lm_robust(fel_arr_bin ~ atl_white_blv + 
                      atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(atl_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')

binmod8 = lm_robust(fel_arr_bin ~ aus_white_blv +
                      aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(aus_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')

binmod9 = lm_robust(fel_arr_bin ~ bos_white_blv + p_race_white +
                      bos_ses_blv + lpop + p_race_black + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                    data = bos_fin, cluster = BG_CODE, se_type = 'stata')

binmod10 = lm_robust(fel_arr_bin ~ chi_white_blv + 
                       chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                       hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                       lpropertycrime + lviolentcrime + phys_bound + com_density, 
                     data = subset(chi_fin, pop > 0), cluster = BG_CODE)

binmod11 = lm_robust(fel_arr_bin ~ mil_white_blv + 
                       mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                       hhi + lmhhi + pown + p_poverty + p_emp_unemployed + 
                       pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                     data = subset(mil_fin, pop > 0), cluster = BG_CODE)

binmod12 = lm_robust(fel_arr_bin ~ sea_white_blv + 
                       sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                       hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                       lpropertycrime + lviolentcrime + phys_bound + com_density, 
                     data = subset(sea_fin, pop > 0), cluster = BG_CODE)

library(modelsummary)
options(modelsummary_format_numeric_latex = 'mathmode')

modelsummary::modelsummary(models = list(negbinm1, negbinm2, negbinm3, 
                                         negbinm4, negbinm5, negbinm6),
                           coef_map = c("atl_white_blv" = "Racial Boundary",
                                        "aus_white_blv" = "Racial Boundary",
                                        "bos_white_blv" = "Racial Boundary",
                                        "chi_white_blv" = "Racial Boundary",
                                        "mil_white_blv" = "Racial Boundary",
                                        "sea_white_blv" = "Racial Boundary",
                                        "p_race_white_blv" = "Boundary (White)",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "p_race_white" = "% White",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "p_poverty" = "% Poverty",
                                        "pown" = "% Homeowner",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "lpropertycrime" = "Log(Property Crime)",
                                        "lviolentcrime" = "Log(Violent Crime)",
                                        "lcrime" = "Log(Crime)",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density"),
                           output = "latex",
                           file = "negbinm_misd.tex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: White/Non-White Racial Geographic Boundaries Are Associated With Misdemeanor Arrests',
                           latex_options = c("scale_down", "tabularx"),
                           file = "negbinm_misd.tex")


modelsummary::modelsummary(models = list(negbinm7, negbinm8, negbinm9, 
                                         negbinm10, negbinm11, negbinm12),
                           coef_map = c("atl_white_blv" = "Racial Boundary",
                                        "aus_white_blv" = "Racial Boundary",
                                        "bos_white_blv" = "Racial Boundary",
                                        "chi_white_blv" = "Racial Boundary",
                                        "mil_white_blv" = "Racial Boundary",
                                        "sea_white_blv" = "Racial Boundary",
                                        "p_race_white_blv" = "Boundary (White)",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "p_race_white" = "% White",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "p_poverty" = "% Poverty",
                                        "pown" = "% Homeowner",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "lpropertycrime" = "Log(Property Crime)",
                                        "lviolentcrime" = "Log(Violent Crime)",
                                        "lcrime" = "Log(Crime)",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density"),
                           output = "latex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: White/Non-White Racial Geographic Boundaries Are Associated With Felony Arrests',
                           file = "negbinm_fel.tex")


texreg(l = list(binmod1, binmod2, binmod3, binmod4, binmod5, binmod6),
       caption = "Binary Outome Model: White/Non-White Racial Geographic Boundaries Are Associated With Misdemeanor Arrests",
       file = "binary_misd.tex",
       custom.coef.map = list('boundary_quart_dummy' = "Racial Boundary",
                              'atl_white_blv' = "Racial Boundary",
                              'aus_white_blv' = "Racial Boundary",
                              'bos_white_blv' = "Racial Boundary",
                              'chi_white_blv' = "Racial Boundary",
                              'mil_white_blv' = "Racial Boundary",
                              'sea_white_blv' = "Racial Boundary",
                              "atl_ses_blv" = "SES Boundary",
                              "aus_ses_blv" = "SES Boundary",
                              "bos_ses_blv" = "SES Boundary",
                              "chi_ses_blv" = "SES Boundary",
                              "mil_ses_blv" = "SES Boundary",
                              "sea_ses_blv" = "SES Boundary",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       custom.model.names = c("ATL", "AUS", "BOS", "CHI", "MIL", "SEA"),
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE)

texreg(l = list(binmod7, binmod8, binmod9, binmod10, binmod11, binmod12),
       caption = "Binary Outome Model: White/Non-White Racial Geographic Boundaries Are Associated With Felony Arrests",
       file = "binary_fel.tex",
       custom.coef.map = list("boundary_quart_dummy" = "Racial Boundary",
                              "atl_ses_blv" = "SES Boundary",
                              "aus_ses_blv" = "SES Boundary",
                              "bos_ses_blv" = "SES Boundary",
                              "chi_ses_blv" = "SES Boundary",
                              "mil_ses_blv" = "SES Boundary",
                              "sea_ses_blv" = "SES Boundary",
                              'atl_white_blv' = "Racial Boundary",
                              'aus_white_blv' = "Racial Boundary",
                              'bos_white_blv' = "Racial Boundary",
                              'chi_white_blv' = "Racial Boundary",
                              'mil_white_blv' = "Racial Boundary",
                              'sea_white_blv' = "Racial Boundary",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       custom.model.names = c("ATL", "AUS", "BOS", "CHI", "MIL", "SEA"),
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE)


#### robustness: using different outcome measures to circumvent problem with log outcome, crime heterogeneity #### 
aus_fin <- aus_fin %>% mutate(boundary_quart = quantile(aus_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))
bos_fin <- bos_fin %>% mutate(boundary_quart = quantile(bos_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))
sea_fin <- sea_fin %>% mutate(boundary_quart = quantile(sea_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

###### moderation between race boundary dummy and total crime; on misdemeanor arrests ###
# neg. binomial models 
nb_crimemod1 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*lcrime +
                                   atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                 data = subset(atl_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod2 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*lcrime +
                                   aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                 data = subset(aus_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod3 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*lcrime + 
                                   bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                                 data = subset(bos_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod4 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*lcrime +
                                   chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                 data = subset(chi_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod5 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*lcrime +
                                   mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                 data = subset(mil_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod6 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*lcrime + 
                                   sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                 data = subset(sea_fin, pop > 0), cluster = ~BG_CODE)

nb_crimemod7 <-fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*lcrime +
                                  atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                data = subset(atl_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod8 <-fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*lcrime +
                                  aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                data = subset(aus_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod9 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*lcrime + 
                                   bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                                 data = subset(bos_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod10 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*lcrime + 
                                    chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                                    hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                                  data = subset(chi_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod11 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*lcrime + 
                                    mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                                    hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                                  data = subset(mil_fin, pop > 0), cluster = ~BG_CODE)
nb_crimemod12 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*lcrime + 
                                    sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                                    hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                  data = subset(sea_fin, pop > 0), cluster = ~BG_CODE)

# binary models
bin_crimemod1 <- lm_robust(mis_arr_bin ~ boundary_quart_dummy*lcrime +
                             atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                           data = subset(atl_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod2 <- lm_robust(mis_arr_bin ~ boundary_quart_dummy*lcrime +
                             aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                           data = subset(aus_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod3 <- lm_robust(mis_arr_bin ~ boundary_quart_dummy*lcrime + 
                             bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                           data = subset(bos_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod4 <- lm_robust(mis_arr_bin ~ boundary_quart_dummy*lcrime + 
                             chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                           data = subset(chi_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod5 <- lm_robust(mis_arr_bin ~ boundary_quart_dummy*lcrime + 
                             mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                           data = subset(mil_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod6 <-lm_robust(mis_arr_bin ~ boundary_quart_dummy*lcrime + 
                            sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                          data = subset(sea_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')

bin_crimemod7 <- lm_robust(fel_arr_bin ~ boundary_quart_dummy*lcrime +
                             atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                           data = subset(atl_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod8 <- lm_robust(fel_arr_bin ~ boundary_quart_dummy*lcrime +
                             aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                           data = subset(aus_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod9 <- lm_robust(fel_arr_bin ~ boundary_quart_dummy*lcrime + 
                             bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                           data = subset(bos_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod10 <- lm_robust(fel_arr_bin ~ boundary_quart_dummy*lcrime + 
                              chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                            data = subset(chi_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod11 <- lm_robust(fel_arr_bin ~ boundary_quart_dummy*lcrime + 
                              mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                            data = subset(mil_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
bin_crimemod12 <- lm_robust(fel_arr_bin ~ boundary_quart_dummy*lcrime + 
                              sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                            data = subset(sea_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')

## save summary tables as .tex files
library(modelsummary)
options(modelsummary_format_numeric_latex = 'mathmode')

modelsummary::modelsummary(models = list(nb_crimemod1, nb_crimemod2, nb_crimemod3, nb_crimemod4, nb_crimemod5, nb_crimemod6),
                           coef_map = c("aus_white_blv" = "Racial Boundary",
                                        "bos_white_blv" = "Racial Boundary",
                                        "sea_white_blv" = "Racial Boundary",
                                        "atl_ses_blv" = "SES Boundary",
                                        "aus_ses_blv" = "SES Boundary",
                                        "bos_ses_blv" = "SES Boundary",
                                        "chi_ses_blv" = "SES Boundary",
                                        "mil_ses_blv" = "SES Boundary",
                                        "sea_ses_blv" = "SES Boundary",
                                        "boundary_quart_dummy" = "Boundary (White)",
                                        "lcrime" = "Log(Crime)",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "p_race_white" = "% White",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "p_poverty" = "% Poverty",
                                        "pown" = "% Homeowner",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density",
                                        "boundary_quart_dummy:lcrime" = "BoundaryXCrime"),
                           output = "latex",
                           file = "nb_crimemod_misd.tex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: Crime Moderation Analysis for White/Non-White Boundaries and Misdemeanor Arrests')


modelsummary::modelsummary(models = list(nb_crimemod7, nb_crimemod8, nb_crimemod9, nb_crimemod10, nb_crimemod11, nb_crimemod12),
                           coef_map = c("aus_white_blv" = "Racial Boundary",
                                        "bos_white_blv" = "Racial Boundary",
                                        "sea_white_blv" = "Racial Boundary",
                                        "atl_white_blv" = "Racial Boundary",
                                        "chi_white_blv" = "Racial Boundary",
                                        "mil_white_blv" = "Racial Boundary",
                                        "atl_ses_blv" = "SES Boundary",
                                        "aus_ses_blv" = "SES Boundary",
                                        "bos_ses_blv" = "SES Boundary",
                                        "chi_ses_blv" = "SES Boundary",
                                        "mil_ses_blv" = "SES Boundary",
                                        "sea_ses_blv" = "SES Boundary",
                                        "boundary_quart_dummy" = "Boundary (White)",
                                        "lcrime" = "Log(Crime)",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "p_race_white" = "% White",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "p_poverty" = "% Poverty",
                                        "pown" = "% Homeowner",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density",
                                        "boundary_quart_dummy:lcrime" = "BoundaryXCrime"),
                           output = "latex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: Crime Moderation Analysis for White/Non-White Boundaries and Felony Arrests',
                           file = "nb_crimemod_fel.tex")


texreg(l = list(bin_crimemod1, bin_crimemod2, bin_crimemod3, bin_crimemod4, bin_crimemod5, bin_crimemod6),
       caption = "Crime Moderation Analysis for White/Non-White Boundaries and Misdemeanor Arrests",
       file = "binary_crimemod_misd.tex",
       custom.coef.map = list("aus_white_blv" = "Racial Boundary",
                              "bos_white_blv" = "Racial Boundary",
                              "sea_white_blv" = "Racial Boundary",
                              "atl_white_blv" = "Racial Boundary",
                              "chi_white_blv" = "Racial Boundary",
                              "mil_white_blv" = "Racial Boundary",
                              "atl_ses_blv" = "SES Boundary",
                              "aus_ses_blv" = "SES Boundary",
                              "bos_ses_blv" = "SES Boundary",
                              "chi_ses_blv" = "SES Boundary",
                              "mil_ses_blv" = "SES Boundary",
                              "sea_ses_blv" = "SES Boundary",
                              "boundary_quart_dummy" = "Boundary (White)",
                              "lcrime" = "Log(Crime)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density",
                              "boundary_quart_dummy:lcrime" = "BoundaryXCrime"),
       custom.model.names = c("ATL", "AUS", "BOS", "CHI", "MIL", "SEA"),
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE)

texreg(l = list(bin_crimemod7, bin_crimemod8, bin_crimemod9, bin_crimemod10, bin_crimemod11, bin_crimemod12),
       caption = "Crime Moderation Analysis for White/Non-White Boundaries and Felony Arrests",
       file = "binary_crimemod_fel.tex",
       custom.coef.map = list("aus_white_blv" = "Racial Boundary",
                              "bos_white_blv" = "Racial Boundary",
                              "sea_white_blv" = "Racial Boundary",
                              "atl_white_blv" = "Racial Boundary",
                              "chi_white_blv" = "Racial Boundary",
                              "mil_white_blv" = "Racial Boundary",
                              "atl_ses_blv" = "SES Boundary",
                              "aus_ses_blv" = "SES Boundary",
                              "bos_ses_blv" = "SES Boundary",
                              "chi_ses_blv" = "SES Boundary",
                              "mil_ses_blv" = "SES Boundary",
                              "sea_ses_blv" = "SES Boundary",
                              "boundary_quart_dummy" = "Boundary (White)",
                              "lcrime" = "Log(Crime)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density",
                              "boundary_quart_dummy:lcrime" = "BoundaryXCrime"),
       custom.model.names = c("ATL", "AUS", "BOS", "CHI", "MIL", "SEA"),
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE)

#### robustness: re-running stops cleaning #### 

## Austin

load("aus_stops_final.RData")

# now run analyses

# log and +1 to variables (pop already logged)
aus_stops_final = aus_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
aus_stops_final = aus_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_stops_final$larrest_sd <- aus_stops_final$larrests - (mean(aus_stops_final$larrests)/sd(aus_stops_final$larrests))
aus_stops_final$lmisdemeanors_sd <- aus_stops_final$lmisdemeanors - (mean(aus_stops_final$lmisdemeanors)/sd(aus_stops_final$lmisdemeanors))
aus_stops_final$lfelonies_sd <- aus_stops_final$lfelonies - (mean(aus_stops_final$lfelonies)/sd(aus_stops_final$lfelonies))
aus_stops_final$lnonviolent_sd <- aus_stops_final$lnonviolent - (mean(aus_stops_final$lnonviolent)/sd(aus_stops_final$lnonviolent))
aus_stops_final$lviolent_sd <- aus_stops_final$lviolent - (mean(aus_stops_final$lviolent)/sd(aus_stops_final$lviolent))
aus_stops_final$lsociety_sd <- aus_stops_final$lsociety - (mean(aus_stops_final$lsociety)/sd(aus_stops_final$lsociety))
aus_stops_final$lperson_sd <- aus_stops_final$lperson - (mean(aus_stops_final$lperson)/sd(aus_stops_final$lperson))
aus_stops_final$lproperty_sd <- aus_stops_final$lproperty - (mean(aus_stops_final$lproperty)/sd(aus_stops_final$lproperty))

# create scaled DVs for stop vars
aus_stops_final$lall_stops_total_sd <- aus_stops_final$lall_stops_total - (mean(aus_stops_final$lall_stops_total)/sd(aus_stops_final$lall_stops_total))
aus_stops_final$lall_stops_black_sd <- aus_stops_final$lall_stops_black - (mean(aus_stops_final$lall_stops_black)/sd(aus_stops_final$lall_stops_black))
aus_stops_final$lall_stops_latino_sd <- aus_stops_final$lall_stops_latino - (mean(aus_stops_final$lall_stops_latino)/sd(aus_stops_final$lall_stops_latino))
aus_stops_final$lall_stops_white_sd <- aus_stops_final$lall_stops_white - (mean(aus_stops_final$lall_stops_white)/sd(aus_stops_final$lall_stops_white))
aus_stops_final$lall_stops_asian_sd <- aus_stops_final$lall_stops_asian - (mean(aus_stops_final$lall_stops_asian)/sd(aus_stops_final$lall_stops_asian))
aus_stops_final$lall_stops_nonwhite_sd <- aus_stops_final$lall_stops_nonwhite - (mean(aus_stops_final$lall_stops_nonwhite)/sd(aus_stops_final$lall_stops_nonwhite))


# create binned boundary measure
aus_stops_final <- aus_stops_final %>% mutate(boundary_quart = quantile(aus_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# race_blvXcrime plots 
## Austin

library(prediction)

mdata_aus <- aus_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy, all_stops_total,
                                              lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                              hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                              com_density, phys_bound) %>% na.omit()


aus1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

## Milwaukee
load("mil_stops_final.RData")


# log and +1 to variables (pop already logged)
mil_stops_final = mil_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
mil_stops_final = mil_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1),
         lped_stops_total = log(ped_stops_total + 1),
         lped_stops_black = log(ped_stops_black + 1),
         lped_stops_latino = log(ped_stops_latino + 1),
         lped_stops_white = log(ped_stops_white + 1),
         lped_stops_asian = log(ped_stops_asian + 1),
         lped_stops_nonwhite = log(ped_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_stops_final$larrest_sd <- mil_stops_final$larrests - (mean(mil_stops_final$larrests)/sd(mil_stops_final$larrests))
mil_stops_final$lmisdemeanors_sd <- mil_stops_final$lmisdemeanors - (mean(mil_stops_final$lmisdemeanors)/sd(mil_stops_final$lmisdemeanors))
mil_stops_final$lfelonies_sd <- mil_stops_final$lfelonies - (mean(mil_stops_final$lfelonies)/sd(mil_stops_final$lfelonies))
mil_stops_final$lnonviolent_sd <- mil_stops_final$lnonviolent - (mean(mil_stops_final$lnonviolent)/sd(mil_stops_final$lnonviolent))
mil_stops_final$lviolent_sd <- mil_stops_final$lviolent - (mean(mil_stops_final$lviolent)/sd(mil_stops_final$lviolent))
mil_stops_final$lsociety_sd <- mil_stops_final$lsociety - (mean(mil_stops_final$lsociety)/sd(mil_stops_final$lsociety))
mil_stops_final$lperson_sd <- mil_stops_final$lperson - (mean(mil_stops_final$lperson)/sd(mil_stops_final$lperson))
mil_stops_final$lproperty_sd <- mil_stops_final$lproperty - (mean(mil_stops_final$lproperty)/sd(mil_stops_final$lproperty))

# create scaled DVs for stop vars
mil_stops_final$lall_stops_total_sd <- mil_stops_final$lall_stops_total - (mean(mil_stops_final$lall_stops_total)/sd(mil_stops_final$lall_stops_total))
mil_stops_final$lall_stops_black_sd <- mil_stops_final$lall_stops_black - (mean(mil_stops_final$lall_stops_black)/sd(mil_stops_final$lall_stops_black))
mil_stops_final$lall_stops_latino_sd <- mil_stops_final$lall_stops_latino - (mean(mil_stops_final$lall_stops_latino)/sd(mil_stops_final$lall_stops_latino))
mil_stops_final$lall_stops_white_sd <- mil_stops_final$lall_stops_white - (mean(mil_stops_final$lall_stops_white)/sd(mil_stops_final$lall_stops_white))
mil_stops_final$lall_stops_asian_sd <- mil_stops_final$lall_stops_asian - (mean(mil_stops_final$lall_stops_asian)/sd(mil_stops_final$lall_stops_asian))
mil_stops_final$lall_stops_nonwhite_sd <- mil_stops_final$lall_stops_nonwhite - (mean(mil_stops_final$lall_stops_nonwhite)/sd(mil_stops_final$lall_stops_nonwhite))


# now create binned racial blv measures
mil_stops_final <- mil_stops_final %>% mutate(boundary_quart = quantile(mil_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))



library(prediction)


mdata_mil <- mil_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy, all_stops_total,
                                              lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                              hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                              com_density, phys_bound) %>% na.omit()


mil1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

## Chicago
load("chi_stops_final.RData")

# now create binned racial blv measures
chi_stops_final <- chi_stops_final %>% mutate(boundary_quart = quantile(chi_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

# now run analyses

# log and +1 to variables (pop already logged)
chi_stops_final = chi_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime + 1),
         lviolentcrime = log(violent_crime + 1))

# now for new stop vars
chi_stops_final = chi_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_stops_final$larrest_sd <- chi_stops_final$larrests - (mean(chi_stops_final$larrests)/sd(chi_stops_final$larrests))
chi_stops_final$lmisdemeanors_sd <- chi_stops_final$lmisdemeanors - (mean(chi_stops_final$lmisdemeanors)/sd(chi_stops_final$lmisdemeanors))
chi_stops_final$lfelonies_sd <- chi_stops_final$lfelonies - (mean(chi_stops_final$lfelonies)/sd(chi_stops_final$lfelonies))
chi_stops_final$lnonviolent_sd <- chi_stops_final$lnonviolent - (mean(chi_stops_final$lnonviolent)/sd(chi_stops_final$lnonviolent))
chi_stops_final$lviolent_sd <- chi_stops_final$lviolent - (mean(chi_stops_final$lviolent)/sd(chi_stops_final$lviolent))
chi_stops_final$lsociety_sd <- chi_stops_final$lsociety - (mean(chi_stops_final$lsociety)/sd(chi_stops_final$lsociety))
chi_stops_final$lperson_sd <- chi_stops_final$lperson - (mean(chi_stops_final$lperson)/sd(chi_stops_final$lperson))
chi_stops_final$lproperty_sd <- chi_stops_final$lproperty - (mean(chi_stops_final$lproperty)/sd(chi_stops_final$lproperty))

# create scaled DVs for stop vars
chi_stops_final$lall_stops_total_sd <- chi_stops_final$lall_stops_total - (mean(chi_stops_final$lall_stops_total)/sd(chi_stops_final$lall_stops_total))
chi_stops_final$lall_stops_black_sd <- chi_stops_final$lall_stops_black - (mean(chi_stops_final$lall_stops_black)/sd(chi_stops_final$lall_stops_black))
chi_stops_final$lall_stops_latino_sd <- chi_stops_final$lall_stops_latino - (mean(chi_stops_final$lall_stops_latino)/sd(chi_stops_final$lall_stops_latino))
chi_stops_final$lall_stops_white_sd <- chi_stops_final$lall_stops_white - (mean(chi_stops_final$lall_stops_white)/sd(chi_stops_final$lall_stops_white))
chi_stops_final$lall_stops_asian_sd <- chi_stops_final$lall_stops_asian - (mean(chi_stops_final$lall_stops_asian)/sd(chi_stops_final$lall_stops_asian))
chi_stops_final$lall_stops_nonwhite_sd <- chi_stops_final$lall_stops_nonwhite - (mean(chi_stops_final$lall_stops_nonwhite)/sd(chi_stops_final$lall_stops_nonwhite))

library(prediction)


mdata_chi <- chi_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy, all_stops_total,
                                              lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                              hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                              com_density, phys_bound) %>% na.omit()


chi1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                 data = mdata_chi, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
          data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)



mdata_aus$stops_bin = ifelse(mdata_aus$all_stops_total > 0, 1, 0)
mdata_mil$stops_bin = ifelse(mdata_mil$all_stops_total > 0, 1, 0)
mdata_chi$stops_bin = ifelse(mdata_chi$all_stops_total > 0, 1, 0)

# neg. binom models
nbinom_stops_crime1 <- fixest::fenegbin(all_stops_total ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                        data = subset(mdata_aus, pop > 0), cluster = ~BG_CODE)

nbinom_stops_crime2 <- fixest::fenegbin(all_stops_total ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                        data = subset(mdata_chi, pop > 0), cluster = ~BG_CODE)

nbinom_stops_crime3 <- fixest::fenegbin(all_stops_total ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density,
                                        data = subset(mdata_mil, pop > 0), cluster = ~BG_CODE)



modelsummary::modelsummary(models = list(nbinom_stops_crime1, nbinom_stops_crime2, nbinom_stops_crime3),
                           coef_map = c("boundary_quart_dummy" = "Boundary (White)",
                                        "lcrime" = "Log(Crime)",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "p_race_white" = "% White",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "pown" = "% Homeowner",
                                        "p_poverty" = "% Poverty",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density",
                                        "boundary_quart_dummy:lcrime" = "Boundary X Crime"),
                           output = "latex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: Influence of Logged Crime on Logged Stops (Standardized), Conditional on Racial Boundary Status.',
                           file = "nb_stops_crimemod.tex")

## binary models 
bin_stops_crime1 <- lm_robust(stops_bin ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                              data = subset(mdata_aus, pop > 0), cluster = BG_CODE)

bin_stops_crime2 <- lm_robust(stops_bin ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                              data = subset(mdata_chi, pop > 0), cluster = BG_CODE)

bin_stops_crime3 <- lm_robust(stops_bin ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                              data = subset(mdata_mil, pop > 0), cluster = BG_CODE)



texreg(l = list(bin_stops_crime1, bin_stops_crime2, bin_stops_crime3),
       caption = "Binary Outome Model: Influence of Logged Crime on Logged Stops (Standardized), Conditional on Racial Boundary Status.",
       file = "bin_stops_crimemod.tex",
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "lcrime" = "Log(Crime)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "pown" = "% Homeowner",
                              "p_poverty" = "% Poverty",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density",
                              "boundary_quart_dummy:lcrime" = "Boundary X Crime"),
       custom.model.names = c("AUS", "CHI", "MIL"),
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE)


#### robustness: using different outcome measures to circumvent problem with log outcome #### 
# Heterogeneous Influence of Racial Boundary on Arrests By White Racial Context. 
form_list = 
  list(
    'lmisdemeanors_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + pown + 
    p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound',
    'lfelonies_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + pown + 
    p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound')

# neg. binom models
# misd
nb_racemod1 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*p_race_white + 
                                  atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(atl_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod2 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*p_race_white + 
                                  aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(aus_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod3 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*p_race_white + 
                                  bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(bos_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod4 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*p_race_white + 
                                  chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(chi_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod5 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*p_race_white + 
                                  mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(mil_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod6 <- fixest::fenegbin(misdemeanor_arrests ~ boundary_quart_dummy*p_race_white + 
                                  sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(sea_fin, pop > 0), cluster = ~BG_CODE)
# fel
nb_racemod7 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*p_race_white + 
                                  atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(atl_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod8 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*p_race_white + 
                                  aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(aus_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod9 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*p_race_white + 
                                  bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = subset(bos_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod10 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*p_race_white + 
                                   chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                 data = subset(chi_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod11 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*p_race_white + 
                                   mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                 data = subset(mil_fin, pop > 0), cluster = ~BG_CODE)
nb_racemod12 <- fixest::fenegbin(felony_arrests ~ boundary_quart_dummy*p_race_white + 
                                   sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                 data = subset(sea_fin, pop > 0), cluster = ~BG_CODE)
# binary models
# misd
binmod1 = lm_robust(mis_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(atl_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod2 = lm_robust(mis_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(aus_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod3 = lm_robust(mis_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lcrime + phys_bound + com_density,
                    data = subset(bos_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod4 = lm_robust(mis_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(chi_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod5 = lm_robust(mis_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(mil_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod6 = lm_robust(mis_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(sea_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
# fel
binmod7 = lm_robust(fel_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(atl_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod8 = lm_robust(fel_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lpropertycrime + lviolentcrime + phys_bound + com_density,
                    data = subset(aus_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod9 = lm_robust(fel_arr_bin ~ boundary_quart_dummy*p_race_white + 
                      bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                      hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                      lcrime + phys_bound + com_density,
                    data = subset(bos_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod10 = lm_robust(fel_arr_bin ~ boundary_quart_dummy*p_race_white + 
                       chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                       hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                       lpropertycrime + lviolentcrime + phys_bound + com_density,
                     data = subset(chi_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod11 = lm_robust(fel_arr_bin ~ boundary_quart_dummy*p_race_white + 
                       mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                       hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                       lpropertycrime + lviolentcrime + phys_bound + com_density,
                     data = subset(mil_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')
binmod12 = lm_robust(fel_arr_bin ~ boundary_quart_dummy*p_race_white + 
                       sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                       hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                       lpropertycrime + lviolentcrime + phys_bound + com_density,
                     data = subset(sea_fin, pop > 0), cluster = BG_CODE, se_type = 'stata')

# print summary tables
modelsummary::modelsummary(models = list(nb_racemod1, nb_racemod2, nb_racemod3, nb_racemod4, nb_racemod5, nb_racemod6),
                           coef_map = c("boundary_quart_dummy" = "Boundary (White)",
                                        "p_race_white" = "% White",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "p_poverty" = "% Poverty",
                                        "pown" = "% Homeowner",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "lpropertycrime" = "Log(Property Crime)",
                                        "lviolentcrime" = "Log(Violent Crime)",
                                        "lcrime" = "Log(Crime)",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density",
                                        "boundary_quart_dummy:p_race_white" = "Boundary X % White"),
                           output = "latex",
                           file = "nb_hetefx_misd.tex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: Influence of Percent White on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status')

# print summary tables
modelsummary::modelsummary(models = list(nb_racemod7, nb_racemod8, nb_racemod9, nb_racemod10, nb_racemod11, nb_racemod12),
                           coef_map = c("boundary_quart_dummy" = "Boundary (White)",
                                        "p_race_white" = "% White",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "p_poverty" = "% Poverty",
                                        "pown" = "% Homeowner",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "lpropertycrime" = "Log(Property Crime)",
                                        "lviolentcrime" = "Log(Violent Crime)",
                                        "lcrime" = "Log(Crime)",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density",
                                        "boundary_quart_dummy:p_race_white" = "Boundary X % White"),
                           output = "latex",
                           file = "nb_hetefx_fel.tex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: Influence of Percent White on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status')



texreg(l = list(binmod1, binmod2, binmod3, binmod4, binmod5, binmod6),
       caption = "Binary Outome Model: Influence of Percent White on Logged Midemeanor Arrests (Standardized), Conditional on Racial Boundary Status",
       file = "bin_hetefx_misd.tex",
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "atl_ses_blv" = "SES Boundary",
                              "aus_ses_blv" = "SES Boundary",
                              "bos_ses_blv" = "SES Boundary",
                              "chi_ses_blv" = "SES Boundary",
                              "mil_ses_blv" = "SES Boundary",
                              "sea_ses_blv" = "SES Boundary",
                              'atl_white_blv' = "Racial Boundary",
                              'aus_white_blv' = "Racial Boundary",
                              'bos_white_blv' = "Racial Boundary",
                              'chi_white_blv' = "Racial Boundary",
                              'mil_white_blv' = "Racial Boundary",
                              'sea_white_blv' = "Racial Boundary",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density",
                              "boundary_quart_dummy:p_race_white" = "Boundary X % White"),
       custom.model.names = c("ATL", "AUS", "BOS", "CHI", "MIL", "SEA"),
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE)

texreg(l = list(binmod7, binmod8, binmod9, binmod10, binmod11, binmod12),
       caption = "Binary Outome Model: Influence of Percent White on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status",
       file = "bin_hetefx_fel.tex",
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "atl_ses_blv" = "SES Boundary",
                              "aus_ses_blv" = "SES Boundary",
                              "bos_ses_blv" = "SES Boundary",
                              "chi_ses_blv" = "SES Boundary",
                              "mil_ses_blv" = "SES Boundary",
                              "sea_ses_blv" = "SES Boundary",
                              "atl_white_blv" = "Boundary (White)",
                              "aus_white_blv" = "Boundary (White)",
                              "bos_white_blv" = "Boundary (White)",
                              "chi_white_blv" = "Boundary (White)",
                              "mil_white_blv" = "Boundary (White)",
                              "sea_white_blv" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density",
                              "boundary_quart_dummy:p_race_white" = "Boundary X % White"),
       custom.model.names = c("ATL", "AUS", "BOS", "CHI", "MIL", "SEA"),
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE)


#### robustness: using different outcome measures to circumvent problem with log outcome #### 
# Heterogeneous Influence of Racial Boundary on Police Stops By Race of Civilian. 

form_list = 
  list(
    'lall_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound',
    'lall_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound'
  )

aus_stops_final$stops_bin_nw = ifelse(aus_stops_final$all_stops_nonwhite > 0, 1, 0)
mil_stops_final$stops_bin_nw = ifelse(mil_stops_final$all_stops_nonwhite > 0, 1, 0)
chi_stops_final$stops_bin_nw = ifelse(chi_stops_final$all_stops_nonwhite > 0, 1, 0)

aus_stops_final$stops_bin_w = ifelse(aus_stops_final$all_stops_white > 0, 1, 0)
mil_stops_final$stops_bin_w = ifelse(mil_stops_final$all_stops_white > 0, 1, 0)
chi_stops_final$stops_bin_w = ifelse(chi_stops_final$all_stops_white > 0, 1, 0)

# neg. binom. models
# stops white
nbinom_stops_race1 <- fixest::fenegbin(all_stops_white ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                         pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + 
                                         com_density + phys_bound,
                                       data = subset(aus_stops_final, pop > 0), cluster = ~BG_CODE)

nbinom_stops_race2 <- fixest::fenegbin(all_stops_white ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                         pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + 
                                         com_density + phys_bound,
                                       data = subset(chi_stops_final, pop > 0), cluster = ~BG_CODE)
nbinom_stops_race3 <- fixest::fenegbin(all_stops_white ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                         pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + 
                                         com_density + phys_bound,
                                       data = subset(mil_stops_final, pop > 0), cluster = ~BG_CODE)


# stops nonwhite
nbinom_stops_race4 <- fixest::fenegbin(all_stops_nonwhite ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                         pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                         com_density + phys_bound,
                                       data = subset(aus_stops_final, pop > 0), cluster = ~BG_CODE)

nbinom_stops_race5 <- fixest::fenegbin(all_stops_nonwhite ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                         pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                         com_density + phys_bound,
                                       data = subset(chi_stops_final, pop > 0), cluster = ~BG_CODE)

nbinom_stops_race6 <- fixest::fenegbin(all_stops_nonwhite ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                         pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                         com_density + phys_bound,
                                       data = subset(mil_stops_final, pop > 0), cluster = ~BG_CODE)



# binary models 
# stops white
bin_stops_race1 <- lm_robust(stops_bin_w ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                               pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                               com_density + phys_bound, 
                             data = subset(aus_stops_final, pop > 0), cluster = BG_CODE)

bin_stops_race2 <- lm_robust(stops_bin_w ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                               pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                               com_density + phys_bound, 
                             data = subset(chi_stops_final, pop > 0), cluster = BG_CODE)

bin_stops_race3 <- lm_robust(stops_bin_w ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                               pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                               com_density + phys_bound, 
                             data = subset(mil_stops_final, pop > 0), cluster = BG_CODE)



# stops nonwhite
bin_stops_race4 <- lm_robust(stops_bin_nw ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                               pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                               com_density + phys_bound, 
                             data = subset(aus_stops_final, pop > 0), cluster = BG_CODE)

bin_stops_race5 <- lm_robust(stops_bin_nw ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                               pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                               com_density + phys_bound, 
                             data = subset(chi_stops_final, pop > 0), cluster = BG_CODE)

bin_stops_race6 <- lm_robust(stops_bin_nw ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                               pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                               com_density + phys_bound, 
                             data = subset(mil_stops_final, pop > 0), cluster = BG_CODE)




# print summary tables
modelsummary::modelsummary(models = list(nbinom_stops_race1, nbinom_stops_race4, nbinom_stops_race2, nbinom_stops_race5,
                                         nbinom_stops_race3, nbinom_stops_race6),
                           coef_map = c("boundary_quart_dummy" = "Boundary (White)",
                                        "p_race_white" = "% White",
                                        "ses_blv" = "Boundary (SES)",
                                        "lpop" = "Log(Population)",
                                        "age_15_35_male" = "Age 15-35 Male",
                                        "hhi" = "Diversity",
                                        "lmhhi" = "Log(MHHI)",
                                        "pown" = "% Homeowner",
                                        "p_poverty" = "% Poverty",
                                        "p_emp_unemployed" = "% Unemployed",
                                        "pcol" = "% College",
                                        "phys_bound" = "Physical Boundary",
                                        "com_density" = "Commercial Density",
                                        "lpropertycrime" = "Log(Property Crime)",
                                        "lviolentcrime" = "Log(Violent Crime)",
                                        "boundary_quart_dummy:p_race_white" = "Boundary * % White"),
                           output = "latex",
                           file = "nb_stops_fow.tex",
                           gof_omit = "rmse|RMSE|R2 Adj|AIC|BIC",
                           stars = TRUE,
                           title = 'Neg. Binomial Model: Influence of Percent White on Police Stops by Race of Civilian, Conditional on Racial Boundary Status')


texreg(l = list(bin_stops_race1, bin_stops_race4, bin_stops_race2, bin_stops_race5, bin_stops_race3,
                bin_stops_race6),
       caption = "Binary Outome Model: Influence of Percent White on Police Stops by Race of Civilian, Conditional on Racial Boundary Status",
       file = "bin_stops_fow.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "pown" = "% Homeowner",
                              "p_poverty" = "% Poverty",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "boundary_quart_dummy:p_race_white" = "Boundary * % White"),
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("Aus: White","AUS: Nonwhite","CHI: White","CHI: Nonwhite",
                              "MIL: White","MIL: Nonwhite"),
       stars = c(0.001, 0.01, 0.05, 0.1))




